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Abstract 



We investigate a method to solve a class of Schrodinger equation eigenvalue problems numerically to very high precision P (from 
thousands to a million of decimals). The memory requirement, and the number of high precision algebraic operations, of the method 
C^) scale essentially linearly with P when only eigenvalues are computed. However, since the algorithms for multiplying high precision 
numbers scale at a rate between P l 6 and P log P log log P, the time requirement of our method increases somewhat faster than P 2 . 
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1. Introduction 



The one-dimensional anharmonic oscillator have been sub- 
C| ject to much investigation since the seminal works by Bender 
Q^and Wu (H [2]| on the behaviour of its perturbation expansion. 
Jh The motivation has often been to extract features and under- 
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standing that can be generalized to more interesting situations, 
like quantum field theories in higher space-time dimensions, or 
to test new approximation methods. 

In this note we report briefly from our work on a class of 
one-dimensional quantum mechanical systems which include 
the one mentioned above, i.e., systems which are modelled by 
Schrodinger equations of the type 



x 2M + 



Jim 



i/,(x) = siftx), (1) 



for some finite (small) M, and real coefficients s and v m . Our 
focus has been on the possible accuracy to which the eigenval- 
ues and eigenfunctions can be found within available computa- 
tional resources (memory and CPU cycles). The algorithm we 
have implemented has modest memory demands; the required 



^ memory scales asymptotically with P like MP, where P is the 
desired precision of eigenvalues or eigenfunctions in decimal 
digits. The number of required algebraic operations (involving 
high-precision numbers) generally also seems to grow asymp- 
totically with P like MP. However, in some cases there is a 
large offset which makes it computationally very expensive to 
obtain P to even a few digits. Further, the time required per 
high-precision algebraic operation (i.e. multiplication or di- 
vision) increases somewhat faster than linearely with P. The 
high-precision numerical library (CLN |5|, built on GMP |6|) 
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we have used has not been parallelized. Thus, our algorithm is 
mostly constrained by wall-clock time. 

Due to space constraints we can in the remainder of this note 
only present examples of our results (sections [2]-[4]l and a brief 
decription of the requirements for obtaining a desired precision 
(section [5]). Solving the differential equations ([T]) numerically 
to very high precision is the most simple and straightforward 
part of our work (section |6j; the analysis of inevitable loss of 
accuracy might be more interesting (section [7]). The behaviour 
of our numerical algorithm is illustrated in sectionJH] 

A more complete description will be given elsewhere Q. 

2. Ground state energy to one million decimals 

As our first proof-of-method we considered the ground state 
of the pure anharmonic oscillator, 



— i//"(x) + x tf/(x) — etff(x), 



(2) 



and computed its ground state energy to 1 000000 + decimals. 
The result, obtained after about 20 days of computing, is 

J, Decimal number 1 
e = 1.060 362 090 484 182 899 647 046 016 692 663 \ 

545 515 208 728 528 977 933 216 245 241 695\ 

943 563 044 344 421 126 896 299 134 671 703\ 

J, Decimal 1 000 
... 304 9 16 644 28 1 633 946 163 324 287 004 26 1 \ 

i Decimal 10 000 (3) 
... 578 044 164777 855 042 412 917 855 188 328\ 

i Decimal 100 000 
... 857 326 052 850 064 563 492 099 229 730 278\ 

J, Decimal 1 000 000 
... 820 139466721 621 064477 821 481 635 914 .. . 
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The full sequence is available on request (for people seeking 
diversion from investigating the numerical patterns of n). 

3. Excited states compared to WKB results 

The result of the previous section may be somewhat unsatis- 
factory to sceptical readers, since there is to our knowledge no 
similar results to compare against. However, it is only slightly 
more challenging to treat excited states. To be fair and square 
we have computed eigenvalue number n = 50000 of equa- 
tion (|2| to 50000 + decimals. For such values of n the WKB 
approximation should be resonably good. Our result is 

J, Decimal number 1 
esoooo = 4024985.730438 698 704313 888 104230563\ 
24 1 82 1 769 405 1 66 607 3 1 3 872 288 \ 
953 655 475 876981 078 813 733 788\ 
i Decimal 50 000 
...545 947155 500441209... (4) 

In comparison, the 12 th order WKB approximation computed 
by Bender et. al. [4] gives 

e*™ 12 =4024985.730438 698 704313 888 104230563\ 
241821769405 166 607 313 872 288\ 
953 657... (5) 

I.e., the relative accuracy of the WKB approximation is 



^WKB-12 
"50000 



- £50 000 



5.163... x 10" 



(6) 



£50000 

From observation of the behaviour of the WKB series we find 
this accuracy to be as expected for the 12 th order approximation 
at this value of n. Some of us plan to return to a more detailed 
analysis of the behaviour of the WKB approximation, which 
may nowadays be extended easily to much higher orders. 

4. Brute force calculation of double-well level-splitting 

Another well analysed situation where we may stress-test our 
method is the calculation of the level splitting between the low- 
est even and odd parity eigenstates of the double-well potential, 



s 2 >]/"(x) + [x 1 - l) 2 \fr{x) = si/r(x), 



(7) 



for small s. The lowest even, el , and odd, e ; , parity states 

.(-) _ J + 
o fc o 



are split by an exponentially small amount Aeo = e q ) ~ e i 



Asymptotically as s — > + , 



Aeo ~ As£ ~ J = 16 



e -4/3s e L(j) 



(8) 



where L(s) = — (j^s -\ ) is given to order s 10 by Zinn-Justi 

0. We have made independent calculations of ei to 30000 + 



We use a different normalization: s = 8g and e = 32gE, where g and E are 
the parameters in (3)- 



digits accuracy for s = 50000 : 

Decimal 28 954 I 



= 0.000 039 999 799 .. . 984 723 697 . 



0.000 039 999 799 .. . 990 905 404 . 



(9) 
(10) 



Equation <|8j agrees with the difference to the expected order, 
Ae^" J - As 



Aen 



1.649... x 10~ 48 ~ 8052 s u . (11) 



The right hand side is of the magnitude expected for the next 
term in L(s). 

We hope the three examples above have convinced the reader 
that it is possible to solve the eigenvalue problems ([T]i to very 
high precision. How precise will of course depend on the pa- 
rameters and which eigenstate we want to investigate. 

5. How to achieve a desired precision 

The eigenvalue condition for equation ([T]i is assumed to be 
that if/(x) — > as x — > +oo. We are unable to impose this exactly 
in our numerical algorithm. However, there is an equivalent 
Robin boundary condition which can be imposed at some finite 
(large) x, 



R(x) 



x M + 



(12) 



We don't know R(x) exactly, but there is for any desired pre- 
cision P a finite value of x such that an approximate R(x) is 
sufficient. The required value of x can be estimated by asymp- 
totic analysis of equation ([T]i as x — > oo (or a WKB approxima- 
tion to include estimates of constant prefactors which cannot be 
found by asymptotic analysis alone). A first estimate is that one 
should choose x so that 

6XP (~ ~s X ^ V(y) ~ S ^ ^ 10 PcS ' W ~ 10 P ' (13) 

if one wants to compute the eigenvalue to P decimals precision. 
Here xq is the largest turning point, and to simplify we have 
assumed that the Robin boundary condition is replaced by a 
Diriclet one, R(x) = oo. 



Equation ( 13 1 is an a priori estimate, which we have tested 



by choosing a very large x to obtain a very accurate eigenvalue 
(so that it may be considered exact), and used this to observe 
the obtainable precision for lower values of x. The obtainable 
precision at a given x is defined as 



P(x) = \g\s{x)-e\ 



(14) 



where s(x) is the eigenvalue found numerically by use of equa- 
tion (|T2j. Figure 1 displays the difference between P e st(^) and 
P(x) in some cases. Note that the difference between the es- 
timated Pest(^) and obtainable P(x) precision varies very little 
with x, and hence can be found numerically from fairly low- 
precision calculations. The difference probably occurs because 



we have neglected a slowly varying prefactor in ( 13 1. 
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Figure 1: The obtainable precision P(x) for various 



eigenvalues e with a Diriclet boundary condition ( 12 1 im- 
posed at x instead of the exact condition. The cases plot- 
ted are (i) for the lowest eigenvalue so of equation Q, 
for which P est (x) ~ 2x 3 /31nl0, (ii) for the 10000 th 
eigenvalue £iqooo> an d (iii) f° r the lowest eigenvalue 
e* w of equation for which P es tO) * 2(x - \f(x + 
2)/(3slnl0). 

It is reassuring that the obtainable precision can be predicted 
to within a few digits: Choosing a too small x leads to the so- 
lution of the wrong eigenvalue problem; choosing a too large 
x leads to a waste of CPU cycles. The obtainable precision 
at fixed x may be improved by some number of decimals by 
using a better Robin boundary condition found by asymptotic 
analysis of equation ([TJ. Such improvement might be useful in 
situations where a few tens of decimals precision is sufficient. 



6. Method for solving the Schrodinger equation 



7. Numerical loss of accuracy 

Note that for the harmonic oscillator the summation < p~5] > 
means computing its ground state if/o(x) = e~ x ~l 2 for large x 
by Taylor expansion. This is certainly not the recommended 
method of computation, due to large cancellations and roundoff 
errors. However, by calculating with sufficient numerical pre- 
cision — which is not prohibitively large — it actually works 
quite well. Further, the location of eigenvalues are determined 
by e- values where tfr(x) = if/(x; s) changes sign very rapidly 
with e. The computation of eigenvalues only is less sensitive to 
cancellations. 

Nevertheless, the effects of roundoff and cancellations must 
be considered. Computing with high-precision floating point 
numbers with D decimals accuracy means that the value of if/(x) 
in equation §15\ is accumulated from numbers with a mantissa 



of (Din 10/ In 2) bits. A contribution of magnitude 10 AD to the 
sum will thus have a round-off error of order 10 AD ~ D . S terms 
of the same magnitude is expected to increase this error by a 
factor VS (with symmetric roundoff), which is an insignificant 
increase. In principle there might also be error amplification 
in the recursion relation, but we have not observed signatures 
of such. Thus, we estimate the numerical accuracy loss to AD 
decimal digits, where 



10 fl 



max { \A m (x)\ ) 



(16) 



The largest term in the sum ( 15 1 is found numerically by moni- 
toring the recursion/summation routine. It may also be a pri- 
ori estimated. We have done the latter in two ways: First, by 
asymptotic analysis of the recursion relations in some simple 
situations (those described in sections [2] and [4]). This analysis 
rapidly becomes complicated. Second, all analytic and numeri- 
cal results found are consistent with the assumption that 



max { \A„,(x)\ } = max ( |^(xe^)| } , 

m id v J 



(17) 



We have postponed the description of our numerical method 
of solving equation ([TJ, due to its (perhaps) disappointingly 
naive simplicity: We make a brute force summation of its Tay- 
lor expansion, 



if,{x) 



2/ 

m=0 



,.2m 



(15) 



where <x = or 1 depending on the parity of the solution. The 
coefficients a m , or more precisely A m (x) = are gen- 

erated recursively from equation ([TJ), starting with Aq(x) = 1 . 
Only the M + 1 last coefficients need to be considered at any 
time while the sum is accumulated; hence the memory require- 
ment is low. 

Since equation (JTJ) has no singular points in the finite plane 
the sum will eventually converge very fast. The number of 
terms N needed in the sum (JT3J can be chosen automatically 
by the recursion/summation routine, but may also be a apri- 
ori estimated. The calculation is done in very-high-precision 
floating-point arithmetic using the CLN C++ library of num- 
bers 0. 



where if/(xe llf ) can be estimated from a WKB-approximation. 
Equation ( 17 i is based on the assumptions that (i) there is al- 



ways a point on the circle xe 1 ^ where cancellations are in- 



significant in the sum (15 1, and (ii) the main contributions to 



the sum come from relatively few terms around the maximum 
termj^] There may be parameter combinations where the first as- 
sumption fails, in which case the equality sign in equation ( fT7| 
should be replaced by > (less helpful for estimations). 
For the example in section[2]we find 



Pest( x ) 



2x i 
3 In 10' 



AD : 



3 In 10 



(18) 



Thus, to compute the ground state to P decimals accuracy we 
must evaluate the wavefunction at x = [QlnlO^P] . If we 



2 Actually, the number of terms contributing to maxj \il/(xe l ' e )\ J is an unim- 
portant correction to AD. For the example in section[2fwe used x = 152, for 
which the estimated maximum is about 10 508386 , and summed less than 10 7 
terms of the Taylor series. Whether the maximum is contributed from one sin- 
gle or almost all terms of the sum makes only a few decimals change in AD. 



3 



want to evaluate the wavefunction to P decimals at this x we 
must choose a numerical precision of D — \P decimals, since 
we will loose AD = \P decimals to roundoff errors. How- 
ever, if we only want to find the eigenvalue e we experience a 
compensating accuracy gain. An uncertainty 6\j/ in the wave- 
function translates to an uncertainty 

<fe = | — ] 8\fi (19) 



in the eigenenergy. In this example (dtfr/ds) as e / 3 as 10^ 2 . 
Hence, the accuracy loss due to roundoff is completely com- 
pensated by the accuracy gain caused by a very large {di/z/ds). 
This implies that our method of computing eigenvalues may 
work well even for ordinary precisions P. There are many cases 
where such complete compensation occur. 
The example in section [4] is different. We find 

^estW^^T^ (x-l) Z (x + 2), 



AD 

dip 
He 



3slnl0 
1 /l 



sin 10 \i X " + 2 X 



(20) 



10 



In this case the accuracy loss is not fully compensated. We note 
that AD as 5/ (6s In 10) (large when s is small) when P is chosen 
small. In this case it will be computationally very expensive to 
obtain results to a few decimals of accuracy by this method. 
However, for very large P the computational cost is similar to 
the example in section[2] The example in section[3]is similar to 
the one section|4] only with more cumbersome expressions. 

8. Numerical observations 
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Figure 2: Some observed behaviour when calculating 
the ground state energy of equation Q to P decimals 
precision: (i) The accuracy loss AD (in decimals) due 
to roundoff error; in this case AD ~ P/2. (ii) The number 
N of terms needed in the Taylor expansion ( 15 I, propor- 
tional to P. (iii) The time T used to evaluate one wave- 
function to required precision, locally we find T ~ P v , 
where v w 2.6 around P = 10000 and v as 2.13 around 
P = 200 000. (Datapoints for P > 300 000 are not gener- 
ated under identical conditions.) 



We monitor many parameters during the numerical compu- 
tations. A sample of those are shown in figure 2, from compu- 
tation of the ground state energy of equation Q. We plot the 
numerically observed accuracy loss AD; this agrees with the a 
priori estimate. This is also observed for all other investigated 
cases. 

Further, the number N of terms required in the sum ( [T5| 
seems to grow linearly with precision P. This is also the case 
for other examples, only with different coefficients of propor- 
tionality. 

The total time to make one evaluation of the wavefunction 
if/(x) also behaves as expected from the number of terms in the 
sum, and the time needed to multiply two high-precision num- 
bers. The drop in computation time near P = 42 000 does not 
seem to be an artifact of variations in the computational en- 
vironment. We believe it is due to a change of multiplication 
algorithm at this precision. 

9. Possible extensions 

The possibilites of extending our method to other systems 
are somewhat limited. It seems straightforward to generalize 
to non-symmetric one-dimensional potentials, to Schrodinger 
equations which have only one regular singular point in the fi- 
nite plane, and to small systems of such equations. We believe 
that systems with two (unseparable) degrees of freedom can be 
constructively approached. 

The evaluation of unnormalized wavefunctions is certainly 
possible, with a time requirement proportional to the number of 
evaluation points. We do not rule out the possibility of comput- 
ing normalized wave functions to high precision. 

Another interesting extension is towards very-high-precision 
computation of Green functions for the same class of models. 
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